import numpy as np


def create(n):
    """
    Create a n-dim Extended Freudenstein and Roth function

    :param n: dimensionality
    :return: x0, f, g, G
    """
    if n < 2:
        raise ValueError("n must be >= 2")
    x0 = np.full(n, -2, dtype=np.float64)
    base = np.arange(0, n - 1, dtype=np.int32)

    def fs(x: np.ndarray):
        assert x.ndim == 1
        assert x.shape[0] == n
        x1 = x[:-1]
        x2 = x[1:]
        f1 = x1 + x2 * ((5 - x2) * x2 - 2) - 13
        f2 = x1 + x2 * ((x2 + 1) * x2 - 14) - 29
        return x2, f1, f2

    def df_dx2(x2):
        d1 = ((5 - x2) * x2 - 2) + x2 * (5 - 2 * x2)
        d2 = ((x2 + 1) * x2 - 14) + x2 * (1 + 2 * x2)
        return d1, d2

    def f(x: np.ndarray):
        _x2, f1, f2 = fs(x)
        return np.sum(f1 ** 2) + np.sum(f2 ** 2)

    def g(x: np.ndarray):
        x2, f1, f2 = fs(x)
        d1, d2 = df_dx2(x2)
        g1 = f1 + f2
        g2 = f1 * d1 + f2 * d2
        grad = np.zeros(n, dtype=np.float64)
        grad[:-1] += g1
        grad[1:] += g2
        return grad * 2

    def h(x: np.ndarray):
        x2, f1, f2 = fs(x)
        d1, d2 = df_dx2(x2)
        h11 = np.full(n - 1, 4, dtype=np.float64)
        h12 = (d1 + d2) * 2
        h22 = (d1 * d1 + f1 * (10 - 6 * x2) + d2 * d2 + f2 * (2 + 6 * x2)) * 2
        he = np.zeros((n, n), dtype=np.float64)
        he[base + 0, base + 0] = h11
        he[base + 1, base + 0] = h12
        he[base + 0, base + 1] = h12
        he[base + 1, base + 1] += h22
        return he

    return x0, f, g, h
